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We examine strain-induced quantized Landau levels in graphene. Specifically, arc-bend strains are found 
to cause nonuniform pseudomagnetic fields. Using an effective Dirac model which describes the low-energy 
physics around the nodal points, we show that several of the key qualitative properties of graphene in a strain- 
induced pseudomagnetic field are different compared to the case of an externally applied physical magnetic 
field. We discuss how using different strain strengths allows us to spatially separate the two components of the 
pseudospinor on the different sublattices of graphene. These results are checked against a tight-binding calcu- 
lation on the graphene honeycomb lattice, which is found to exhibit all the features described. Furthermore, we 
find that introducing a Hubbard repulsion on the mean-field level induces a measurable polarization difference 
between the A and the B sublattices, which provides an independent experimental test of the theory presented 
here. 



In the absence of other instabilities, electronic materials in a homogeneous external magnetic field are known to form quantized 
Hall states whose wave functions are localized at the edges of the system. In the usual Hall bar configuration, these quasi-one- 
dimensional states, typically visualized with semicircular classical trajectories, carry opposite currents along the two edges and 
are characterized by the bulks nontrivial Chern number. 

Despite arising from the interaction of noninteracting electrons with a magnetic field, Landau levels can be used to account 
for many of the features of the integer QHE. For example, the large degeneracy of each Landau level (at large magnetic fields) 
accounts for the plateaus in resistivity measurements. For a uniform magnetic field, the existence of propagating states only at 
the edges of the sample can be understood heuristically in terms of the semiclassical trajectories of the electrons near the edges 
of the sample. Therefore, a characterization of the Landau levels becomes a valuable tool to predict the QHE. Furthermore, the 
introduction of a gradient in the applied magnetic field can in principle introduce propagating bulk states, 1 as illustrated in Fig.[T] 

Introducing a strain in graphene quantizes the energy levels'^ and these quantized energy levels mimic the organization of 
the Landau levels. For this reason, the language "pseudo"-magnetic field is used to describe the effect of the strain. In light 
of the importance of Landau levels to the QHE, strain has been proposed as a useful means to study a strain-induced QHE in 
graphene as well as controlling the properties of the electrons in graphene.^ 

In this work, we explore the possibility of using strain in graphene to study the Landau level structure in the presence of 
an effective nonuniform magnetic field. We do this using the Dirac equation coupled to a global U(l) gauge field, which has 
been shown to be an appropriate effective theory near the K and K' point of graphene. 5 - We focus on the case of a linearly 
varying magnetic field, starting from a nonzero value on one side of the quantum Hall bar and decreasing in magnitude to on 
the opposite side. We first show that the only propagating states in the lowest Landau level (LLL) are edge states, and the first 
instance of chiral modes propagating in the bulk occurs in the first excited Landau level. This is in contrast to the nonrelativistic 
case^ where the LLL contains propagating edge and propagating bulk states. We then demonstrate that for a strain-induced 
pseudomagnetic field, only a single zigzag edge of the graphene ribbon supports propagating edge modes, unlike the case of an 
externally applied physical magnetic field that generates oppositely propagating modes at the two zigzag edges. Furthermore, 
since the two components of the Dirac spinor physically correspond to the wave function on the A and B sublattices of graphene, 
we find that changing the magnetic field gradient allows us to observe a physical separation between the two spinor components 
(and hence the wave functions on the A and B sublattices). These results are then confirmed by a tight-binding calculation, 
where we also show how Hubbard repulsion gives rise to opposite magnetic polarizations in the A and B components of the 
spinor, thus making their experimental observation using magnetic tunneling tips feasible. 



We begin by examining the properties of graphene subject to an externally applied physical magnetic field with a finite field 
gradient, both as a review and to emphasize the differences observed when studying the behavior of graphene in strain-induced 



I. INTRODUCTION 



II. LOW-ENERGY PHYSICS IN TERMS OF DIRAC EQUATIONS 



pseudomagnetic fields. The low-energy physics near the nodal K 
described using a set of Dirac equations, 7 : 
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FIG. 1. Semiclassical trajectories of the propagating states in a two-dimensional electron gas in a gradient magnetic field. Besides the 
semicircular edge states, also present in the case of a constant magnetic field, there is a bulk propagating state. 

where vp = -f a a /h, j sa 3eV, a = 1.42A, and e < is the electron charge. Here we focus on a graphene ribbon geometry 
extended in the x direction, such that one can use a plane wave ansatz in that direction and decouple the equations for wave 
functions on the two sublattices at K with vf" = (<J>a, $b) exp(ik x x). This yields a system, 

{-d 2 y + (-k x + eA x f + edyA^ <S> A = e 2 $ A , (2) 

-d 2 y + {-k x + eA x ) 2 - ed y A^j $ B = e 2 $ B , (3) 

and an analogous set of equations at K' with ^' = ($>' A , &' B ) exp(ik x x), where we have chosen the Landau gauge such that 
only A x is non-zero. 8 There is a unique boundary condition for the pseudospinor components, which is determined by the type 
of edge termination in the graphene ribbonP 

The most frequently studied termination types for graphe ne ribbons are armchair and zigzag terminations, although there 
have been many studies on different types (see e.g. refs! 10 ^ 13 !). At zero magnetic field, armchair terminations do not lead to edge 
states, 14 while when a magnetic field is applied, edge states appear at all possible terminations, with the strongest edge states 
at zigzag terminations^. For this reason, we only consider zigzag terminations at y = and y = L, leading to the boundary 
conditions 

$A(y = 0) = <£' A (y = 0) =0, $ B (y = L) = <!>' B (y = L)=Q. (4) 

The differential equations for and Q' B can be solved numerically, subject to their respective boundary conditions. $b and 
Q' A are then given by 

<t> B = ~ kx + eAx ~ dy <t> A , = kx ~ eAx + °y . ( 5) 

e e 
Assuming a linearly varying magnetic field, we parametrize the gauge field profile as 

A x = -^Hy(l-L). (6) 



v F \e\ V2 

Furthermore, it is convenient to define dimensionless coordinates and fields given by 

1 ~ To ~ 

y = a y, L = a L , k x = —k x , e = 7 e , A x = — — A x . (7) 

a \e\VF 

In terms of these dimensionless units, the equations of motion are invariant under the scalings 

y H 

y ->• — , k x -» ak x , e -> as , H -> — , (8) 
a a 

so the scale invariant quantity one should consider is the ratio of H/L. 

In Fig^a), we show the energy spectra for the two LLLs. Similarly to the case of a uniform magnetic field, one observes the 
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coexistence of highly dispersive and practically nondispersive regions in momentum space. Note that the spatial variance of the 
nonuniform magnetic field studied here does not affect this momentum direction. For the LLL(n = 0) around the K and K' 
points, the dispersionful edge states mimic the behavior observed in the two-dimensional electron gas (2DEG), except that in 
graphene the K and K' points provide a single edge state each. This is similarly true for the case of a uniform magnetic field, 
and it is the natural generalization of the single-edge result where only one of the nodal points contributes an edge stated The 
corresponding spatial profiles of the wave functions are shown in Fig. [2^b^. Around the K point, the exponentially decaying 
edge state is at y = and resides exclusively on the B sublattice, whereas around the K' point, the exponentially decaying edge 
state is at y = L and resides exclusively on the A sublattice (which is also the case for a uniform magnetic field). Furthermore, 
the spatial extent of the edges states at each end of the ribbon is different; the edge state at y = L has a larger spatial extent, 
which fits nicely with the semi-classical picture of having a larger half orbit at the edge with a lower magnitude magnetic field 
(which is not the case for a uniform magnetic field). Finally, a propagating bulk mode is observed in the next higher energy level 
(n = 1), also shown in Fig. J^b). 




FIG. 2. (a) Energy spectrum around the K and K 1 points for H = 1 and L = 10 [denned in Eq. |7|] in the presence of a linearly varying 
external magnetic field. The solid (dashed) red curves correspond to the lowest Landau level (n — 0) around the K (K') point, and the solid 
(dashed) blue curves correspond to the n — 1 level around the K (K') point, (b) Magnitudes of the wave functions of the n = edge states 
and n — 1 bulk state. The blue curves around y — are the LLL evaluated at k x — —6.6297 and the red curves around y = L are the LLL 
evaluated at k x = —48.9882 corresponding to a dimensionless energy [Eq. of e 2 w 0.22. The curves around y = 4 are n = 1 bulk states 
evaluated at k x — —30 around the K point (solid (blue) curves) and K' point [dashed (red) curves] corresponding to a dimensionless energy 
[Eq. {8}] of e 2 « 12.61. 

We now turn our attention to the effects of a strain-induced pseudomagnetic field in graphene. The induced pseudogauge field 
alters the equations of motion at the K and K' points to 

(-a x ,-a y )-(p-A)y = —V, (a X) -a y )- (p + X) & = — (9) 
v / vp \ / vp 

We assume a similar form for the gauge field profile as in the physical magnetic field: 

A x = -^Hy(l-L) . (10) 
vp \2 J 

We show in appendix [B] that for a weak strain as depicted in Fig. [5] to a good approximation we can achieve a linear gradient 
pseudomagnetic field as captured by the above gauge field. In order to solve for the energy spectrum, we can proceed in a similar 
fashion as for the physical magnetic field. Results are presented in Fig. |3ja). Since the pseudomagnetic field does not break 
time reversal symmetry, the K and K' spectra are symmetric under sign reversal of k x in this case. The spectrum is separated 
into discrete levels, and the LLL again has a dispersionless regime for a broad range of k x centered around k x = 0, which 
increases in size as the strain is increased. The n = 1 level again captures the propagating bulk mode, although the two K points 
provide oppositely propagating bulk modes. Focusing our attention on the LLL edge states, the propagating edge states occur 
for large \k x \ and have opposite group velocities. In contrast to the case of an externally applied physical magnetic field, we find 
that the propagating edge states in the pseudomagnetic field exist only at the y = L boundary and solely on the A sublattice 
[shown in Fig. |3jb)]. This is a consequence of the fact that the pseudomagnetic field does not break time -reversal symmetry; 
since the spectrum allows for only one edge state (propagating in opposite directions) from each nodal point, time-reversal 
symmetry requires that both modes exist at the same edge. In order to observe the edge state residing around y = (and on the 
B sublattice), one has to move deep into the dispersionless part of the spectrum. As one moves into the dispersionless part of 
the spectrum, i.e. towards lower values of \k x \, the wave function moves from the y = L edge towards the y = edge. In 
this process, the $a wave function grows sharper at the y = L edge. 

The behavior of $s as we move into the dispersionless regime results in a noticeable effect on the local density of state 
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FIG. 3. (a) Energy spectrum around the K and K' points for H — 1, L = 10 [defined in Eq. q7J»] in the presence of a strain induced 
pseudomagnetic field. The solid (dashed) red curves correspond to the lowest Landau level (n = 0) around the K (K') point, and the solid 
(dashed) blue curves correspond to the n = 1 level around the K (K 1 ) point, (b) Magnitude of the wave functions of the LLL edge states 
and n — 1 bulk state. The curves around y — are the LLL evaluated at k x = 48.98707 (blue curves) and k x = —48.98707 (red curves) 
corresponding to a dimensionless energy of e 2 « 0.22. The curves around y = 4 are n = 1 bulk states evaluated at k x = — 30 around the K 
point [solid (blue) curves] and K' point [dashed (red) curves] corresponding to a dimensionless energy [Eq. {8}] ofe 2 « 12.61. 

(LDOS) 17 as the strain is changed. To illustrate this point, we consider all states in an energy window e £ [0, e]. We know 
already that for e m 0, the $a wave function is strongly peaked around y = L (Fig. |4|a)), and for sufficiently small k x , the 
$s wave function is strongly peaked around y = . These states give rise to an enhanced LDOS at the edges. Let us now 
examine the contributions to the LDOS due to the states above zero energy, i.e. the dispersionful states. As shown in Fig. |4jb), 
when the field gradient is increased, the peak of this LDOS shifts towards the y = L edge and grows in strength. Therefore, for 
sufficiently strong strains this feature manifests itself as a distinct peak away from the A sublattice edge state. 




FIG. 4. Contributions to the local density of states from the K and K' point states with e 2 € (0, 0.25] for four different gradient pseudomag- 
netic fields, (a) LDOS on the A sublattice. (b) LDOS on the B sublattice. 



III. LATTICE MODEL CALCULATIONS 

In order to test whether the effective low-energy Dirac model accurately describes the physics of graphene under strain, we 
now turn to the Hubbard model on the honeycomb lattice, which has been extensively used to study the magnetic response of 
graphene (see e.g. refsP^iD), to scrutinize the above results. Specifically, we model the electronic properties of the graphene 
sheet using the Hamiltonian 

N 

H= X! *Pf"r( C ia C J> +h.C.) + U^2 n i<r n i-<r ' C 11 ) 

where the first term is the single-orbital tight-binding model for graphene (with only nearest-neighbor hopping), and the second 
term is the on-site Coulomb repulsion between electrons of opposite spin. Cj CT and cJ CT are electron annihilation and creation 
operators respectively, tp p7T is the hopping integral, and is the number operator given by c\ a Ci a . The coefficient U is chosen 
to be either zero (i.e. a pure tight binding model) or U = l.2t ppn , which is below the critical value of U cr /t pp7I — 2.2 at which 
there is a quantum phase transition to an antiferromagnetic insulating phas d 21 * 2 ^ . For these values of the Hubbard repulsion, the 
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system is in a semi -metallic phase with a conical dispersion, as is experimentally observed in graphene 23 . In the case of a finite 
Hubbard repulsion, the interaction term is simplified using a self-consistent mean field approximation 24 . More details of this 
procedure can be found in appendix [A] 

Applying a strain to the graphene sheet displaces the carbon atoms from their equilibrium positions, and the resulting change 
in bond lengths can be calculated using^l 

5a = ~— Ra ■ (Ui - Ui+a) , 0=1,2,3, (12) 

a 

where ui is the displacement from equilibrium of the i-th A atom at position (x, y) and u i+a is the displacement from equilibrium 
of the neighboring three B atoms. In turn, the change in the bond length modifies the hopping integral tp pv tcPSl 

i^ = - 70 e- 3 - 37 ^. (13) 

In order to generate an approximately linear gradient pseudomagnetic field, as discussed in the previous section, we consider 
the very particular inhomogeneous arc -bend strain 27 illustrated in Fig. |5](here the setup is rotated by 90°). More details can be 
found in appendix [B] The numerical calculations discussed below are performed in real space on a finite 4800-site honeycomb 
lattice. When the system is finite in the y direction, we find no significant spatial variation in that direction for the lowest energy 
states in which we are interested, and therefore, to compute the spectrum as a function of momentum we use periodic boundary 
conditions in the y direction. 




FIG. 5. Illustration of an inhomogeneous arc-bend strain, leading to an approximately constant gradient in the induced pseudomagnetic field. 
An originally rectangular graphene ribbon of length L and width W is distorted into a fan-shaped segment of a circular shell with inner radius 
R and outer radius R + L. 

In Fig. [6] the energy spectra of the strained and unstrained graphene lattice are compared. In the absence of an applied strain, 
the spectrum exhibits a dispersionless band between the two nodal points and a subsequent rise in the dispersion near the nodal 
points that give rise to edge states at zigzag edges 1 - . Upon application of the strain, the dispersionless part of the spectrum 
grows and the behavior around the K and K' points matches the qualitative predictions of the effective low-energy Dirac model 
discussed above. 

In Fig. [7] we show the LDOS as a function of position and energy. It is observed that for small strains, the LDOS remains 
peaked at the edges, but as the strain is increased, a bulk LDOS begins to emerge. On closer analysis (Fig. |7|c-d), we can see 
that this bulk LDOS increases and moves towards the stretched edge as the strain is increased. Furthermore, the LDOS at the 
stretched edge grows in intensity while the LDOS at the smaller edge decreases. These observations are consistent with the 
low-energy continuum model. 

Finally, to show that indeed a bulk peak in the LDOS exists on the B sublattice, we consider the case of a finite Hubbard 
repulsion. We find that turning on U/t > does not change the features of the results so far, except that the system now exhibits 
a local polarization. We show the local polarization in Fig. [8] Here it is observed that the stretched outer edge is associated 
with a negative polarization, whereas the smaller inner edge and the bulk mode are associated with a positive polarization, 
demonstrating that indeed the smaller edge and bulk mode exist on the same sublattice. 

As observed in the effective low-energy continuum model, the bulk mode and the stretched edge mode are actually part of the 
same dispersionful solution of the Dirac equation, and it is therefore very interesting to see such a distinct separation between 
the two. 



6 



(a) 



K' 


/\ K 


j 






2 






/ 1 


1.0 -y 


k -0.6 -J 


\ -0.2 






X 1 






X. -2 










FIG. 6. The two lowest-lying energy levels calculated using the tight-binding model (U = 0). (a) Energy spectrum of the tight-binding model 
on the graphene lattice without a strain, (b) Energy spectrum of the tight-binding model on the graphene lattice in the presence of an arc-bend 
strain with R = 5W. 
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FIG. 7. The local density of states (LDOS) as a function of energy in a graphene sheet with arc-bend strain, calculated using the tight-binding 
model (U = 0). (a ,b) Darker shades correspond to a higher intensity LDOS. (c, d), LDOS near the edges in the tight-binding mode of a 
graphene ribbon is shown. The blue curve is for an arc -bend strain of R = 5.1W, the orange curve corresponds to R = 3.6W, and the red 
curve corresponds to R — 3W.(a) Arc-bend strain with R = 5W. (b) Arc -bend strain with R — 3W. (c) LDOS near smaller inner edge 
integrated over a small energy window around e — 0, (d) LDOS near stretched outer edge integrated over a small energy window around 
e = 0. 



IV. CONCLUSIONS 



We have discussed how strain-induced non-uniform pseudomagnetic fields in graphene lead to propagating quantum Hall 
modes with non-trivial properties. In particular, we considered the case of an approximately linearly varying pseudomagnetic 
field, caused by an arc -bend strain. In contrast to the case of an externally applied magnetic field, the propagating pseudomag- 
netic field edge modes are predominantly localized on the stretched outer zigzag edge of the graphene ribbon. In addition, we 
observe that we can spatially resolve the A and B components of the pseudospinor for this propagating mode with increasing 
arc -bend strain, whereby the B-component appears as a peak in the bulk (as opposed to the edge as in the A case). We also 
find that a finite Hubbard repulsion leads to a polarization of the propagating modes, which makes them accessible to spatially 
resolved magnetic measurements. 
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(a) (b) 




FIG. 8. Local magnetization in a graphene sheet calculated using the Hubbard model. Red symbols denote negative z polarization, whereas 
blue symbols denote positive z polarization, (a) Unstrained case, (b) Strained case. 



In order to evaluate the feasibility of observing this phenomenology experimentally, let us now discuss some estimates. In 
the arc -bend scenario we considered, the strain percentage is given by W/2R, and the maximum strain occurs at the outer 
edge. For example, if R = 3W, the maximum strain is about 17%, and for R = 5W, it would be 10%. Previous theoretical 
papers have predicted that strains in graphene up to 20% can be achieved 28 , whereas so far experimental papers have reported 
strains up to 12%l 2 ^l Therefore, the regime of strains that were explored in our model calculations are of the same order of 
magnitude of what can be realized experimentally. Furthermore, we expect that a similar phenomenology, i.e. strain induced 
polarized edge and bulk modes, should be observable for other strain geometries which may be easier to realize, such as local 
deformation of the lattice. In this context, we would be particularly interested in engineering strain profiles that allow for both 
positive and negative induced pseudomagnetic fields, which would permit us to study the snake mode scenario which had been 
proposecPSEH. Furthermore, the strain-induced vector potential we have studied here accounts for the hopping perturbation only. 
There are corrections to the vector potential arising from deformations of the Brillouin zone that are not the same at the K and 
K' points^, and it would be interesting to study how the physics discussed here is corrected by these terms. 
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Appendix A: Self-Consistent Mean Field Approximation 

Within mean field theory, the repulsion term in the Hubbard model (JTTJ can be approximated by 

N N 

U ^2 n ia ni-cr = U^(- (n la ) (rii_ CT ) + (n icr ) raj_ CT + n. lrJ {m-„)) , (Al) 

i i 

where the mean field (n^) is computed self-consistently from 

(nur) = / dEg zr7 (E)f(E - Ej) . (A2) 

Here gi„{E) = Y^j ^*i{Ej)^i{Ej)5(E — Ej) is the local electronic density of states, Ej is the j'-th energy eigenvalue, and 
f(E — Ef) is the Fermi function. The self-consistent solution provides the local densities of states and the spin densities M* on 
each lattice site, given by 

Mi = ((n m ) - (ni_ CT ))/2 . (A3) 



Appendix B: Strain-Induced Pseudomagnetic Field 

Deforming each lattice point in graphene from position (x, y) to position (x + u x ,y + u y ) leads to an arc -bend strairP^I 

( W\ L 

x + u x = I R + x - — J cos 9(y) - R + — , 

y + u y = (r + x- sm0(y) , (Bl) 
where 9(y) — 2y arcsin(W / /2i?) /W, and R is the radius of the inner side. This induces a gauge field given bjE^l 

A x = -2c—d y u x , A y = c— (d x u x - d v u y ) , (B2) 
where f3 = — d ln(7 ) / d ln(a ) ~ 2. The resulting pseudomagnetic field in the z direction is then given by 

B * = ~w^ cos9{v) arcsin \m ) V ~ w [ R + x ~ 2 ) arcsin ^2H ) ) (B3) 

In the limit of weak strain (W/R —> 0), the pseudomagnetic field is a constant, 

lim B z = -^-, (B4) 
w/R-yO Ra a 

whereas for finite but weak strains the induced pseudomagnetic field is approximately linear in x with maximum magnitude at 
x = and minimum magnitude at x = L. 
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